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Abstract — Astrometric  and  photometric  data  fusion  for  the 
purposes  of  simultaneous  position,  velocity,  attitude,  and  angular 
rate  estimation  has  been  demonstrated  in  the  past.  This  state 
estimation  is  extended  to  include  the  various  surface  parameters 
associated  with  the  bidirectional  reflectance  distribution  function 
(BRDF).  Additionally,  a  physically  consistent  BRDF  and 
radiation  pressure  model  is  utilized  thus  enabling  an  accurate 
physical  link  between  the  observed  photometric  brightness  and 
the  attitudinal  dynamics  and  ultimately  the  orbital  dynamics.  An 
example  scenario  is  then  presented  where  the  model  is  an 
uncontrolled  High  Area  to  Mass  Ratio  (HAMR)  object  in 
geosynchronous  Earth  orbit  and  the  position,  velocity,  attitude, 
angular  rates,  and  surface  parameters  are  estimated 
simultaneously 
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I.  Introduction 

Wetterer  and  Jah  [1]  first  demonstrated  how  brightness 
(photometric  flux  intensity)  measurements  can  be  used  to 
estimate  the  attitude  and  angular  rates  of  a  space  object  (SO). 
Linares  et  al.  [2]  recently  demonstrated  a  mechanism  that 
fuses  both  angles  (line-of-sight)  and  brightness  measurements 
for  the  purpose  of  orbit,  attitude,  and  shape 
determination.  Whereas  angles  measurements  are  direct 
observations  of  the  SO’s  orbital  position,  the  brightness 
measurement  is  dependent  on  orbital  position  and  the  SO’s 
shape,  surface  parameters  and  attitude,  and  thus  provides  an 
indirect  observation  of  these  other  attributes.  The  physical 
correlation  between  the  SO’s  shape/attitude  and  its  orbital 
position  are  caused  by  the  various  non- gravitational  forces  and 
torques,  such  as  radiation  pressures  that  produce  a  linear  and 
angular  acceleration  on  the  SO.  These  radiation  pressures 
must  be  consistent  with  the  surface  bidirectional  reflectance 


distribution  functions  (BRDFs),  as  shown  by  Wetterer  et  al. 

[3]. 

The  work  presented  here  expands  the  parameters  that  are 
included  in  the  state  to  include  those  associated  with  the  space 
object’s  surface,  namely,  the  various  BRDF  parameters.  First, 
the  parameters  included  in  the  SO’s  “augmented”  state  are 
presented  and  the  state  function  detailing  the  dynamics  for  each 
of  these  parameters  in  how  they  are  propagated  forward  in  time 
is  defined.  It  is  important  to  note  that  through  radiation 
pressure,  the  BRDF  parameters  and  other  potential  parameters 
in  the  augmented  state  influence  the  dynamics  that  affect  both 
the  attitude  and  orbit  and  thus  influence  the  evolution  of  the 
system  from  point  to  point.  As  with  the  “classical”  state 
parameters,  these  new  parameters  are  intrinsic  and  unique 
properties  of  the  system.  Next,  the  unscented  Kalman  Filter 
(UKF)  is  briefly  reviewed.  Finally,  test  cases  using  an 
example  scenario  are  examined  in  detail  to  demonstrate 
simultaneous  estimation,  effect  of  model  mismatch,  and 
information  dilution. 

II.  Building  the  Augmented  State 

An  Earth  orbiting  SO’s  attitude,  angular  rates,  position, 
velocity,  size,  shape,  mass  and  surface  characteristics  are  all 
needed  for  high  fidelity  orbit  propagation  and  in  calculating 
associated  measurements  that  are  remotely  observable,  such  as 
angles  and  brightness.  All  of  these  parameters  make  up  the 
SO’s  “augmented”  state.  In  this  paper,  we  will  include  the 
SO’s  attitude,  angular  rates,  position,  velocity,  and  surface 
characteristics  in  the  augmented  state. 

In  this  paper,  the  quaternion,  which  is  based  on  the  Euler 
angle/axis  parameterization  and  contains  four  values,  is  used  to 
specify  the  SO’s  attitude.  The  quaternion  is  defined 
as  q  =  U  With  s  =  [9l  q2  <jr3  ]  =  e  sin(u/2)  and 
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qA  =  cos(l>/2),  where  e  and  v  are  the  Euler  axis  of  rotation  and 
rotation  angle,  respectively.  In  addition,  each  component’s 
angular  rate  is  denoted  by  co  =  [ cox  coy  coz ]  as  defined  in  the 
SO’s  body- fixed  frame.  A  single  position  and  velocity 
corresponding  to  the  SO’s  center  of  mass,  denoted  by  /  =  \x  y 
z]1  and  v1  =  [vx  vy  vz]7  respectively,  are  used,  where  the 
superscript  /  indicates  the  inertial  frame.  This  is  the  classic 
6DOF  representation  of  the  SO’s  orbit. 

Each  material  that  makes  up  the  SO  could  reflect  light 
differently.  Thus,  a  SO  might  have  many  parameters  that 
specify  its  surface  properties.  The  function  that  defines  how 
light  is  reflected  from  an  opaque  surface  with  a  given  surface 
normal  direction  (  N ),  illumination  direction  (  L  with  angles  6t 
and  </>i  from  N ),  and  observer  direction  (  V  with  angles  0r  and 
(j)r  from  N  )  as  shown  in  Fig.  1  is  called  the  bidirectional 
reflectance  distribution  function  (BRDF). 


A 


Fig.  1 .  The  geometry  of  reflection 

The  BRDF  is  given  by 
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where  dLr  is  the  reflected  radiance  in  Wm"2sr-1  and  dEt  is  the 
irradiance  in  Wm"2.  The  bisector  vector  between  the 
illumination  source  and  the  observer  is 


each  to  the  total  ( d  and  s  respectively  where  d  +  s  =  1).  These 
bidirectional  reflectances  are  calculated  differently  for  the 
various  models.  In  this  paper  we  will  use  the  Ashikhmin- 
Shirley  BRDF  [4],  also  known  as  the  Anisotropic  Phong 
BRDF,  where  the  diffuse  and  specular  bidirectional 
reflectances  are  calculated  using 
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where  (4)  is  a  non-Fambertian  diffuse  BRDF,  and  the  Fresnel 
reflectance  (F)  in  (5)  is  given  by  Schlick’s  approximation  [5] 


F  =  F0+fX-FoJi-V.Hj  (6) 

In  addition  to  d ,  p  and  F(h  the  Ashikhmin- Shirley  BRDF  has 
two  exponential  factors  ( nw  nv )  that  define  the  anisotropic 
reflectance  properties  of  each  surface.  Without  loss  of 
functionality,  the  diffuse  reflectance  and  the  specular 
reflectance  at  normal  incidence  can  be  set  equal  to  each  other 
( p  =  F0)  and  the  difference  between  the  diffuse  and  specular 
reflectances  displayed  in  the  diffuse  fraction  parameter,  d. 
Additionally,  for  the  sake  of  simplicity,  in  this  paper  the  two 
exponential  factors  are  set  equal  to  each  other  as  well  ( nu  =  nv 
=  n).  Thus,  there  are  three  unique  surface  parameters  per 
surface  («,  p,  d).  These  surface  parameters  have  constraints  ( n 
>0,  0</)<l,0<J<l).  To  account  for  the  constraints  within 
the  filter,  unconstrained  proxy  values  are  used  in  the  state 
vector  and  estimation  filter  and  these  proxy  values  are 
converted  back  to  the  surface  parameter  value  when  needed. 
The  conversion  equations  to  the  proxy  value  and  from  the 
proxy  value  for  each  of  the  surface  parameters  are 


Pi=  ln(«)>  n  =  exp(p, ) 


(7) 


h  =  (l+v)/|l+v 


(2) 


with  angles  a  and  /?  from  N  and  is  used  in  many  analytic 
BRDF  models. 

There  are  many  different  reflectance  models  that  could  be 
used,  but  all  can  be  expressed  in  a  common  nomenclature  with 
the  general  BRDF  calculated  by 


fr  =  (dRj  +  sRs )  (3) 

which  depends  on  the  diffuse  bidirectional  reflectance  ( Rd )  and 
the  specular  bidirectional  reflectance  ( Rs )  and  the  fraction  of 


p 

i -p 


yC-  =  t-(tanh(/?2)+1) 


(8) 


Pi  =  \ ln( jry]’  d=\ (tanh0>3 )+  0  (9) 

The  shape  model  will  be  built  as  the  sum  of  facets  where 
each  facet  has  a  position  in  the  body  frame  and  is  specified  by 
a  particular  area  and  normal  vector.  The  SO’s  brightness  is 
also  calculated  by  summing  the  contribution  to  the  brightness 
by  each  facet  using  the  BRDF  in  the  equation 


mobjec,=mSun  -2-51og10 


(16) 


Z“4(/J,(n,-l)(n,.-v) 


(10) 


where  At  is  the  area  of  the  zth  facet,  r  is  the  distance  between  the 
SO  and  observer,  and  mSun  is  the  apparent  magnitude  of  the 
illumination  source  (in  this  case  the  Sun). 

So,  in  summary,  the  augmented  state  is 


x  =  [q  co  r7  v7  p]r  (11) 

where  x  is  a  1x16  vector  of  real  numbers  and  p  =  \p\  P2P3]  is 
a  vector  containing  the  surface  parameter  proxy  values. 

III.  Parameter  Propogation 

The  Newtonian  two-body  gravitational  equations  of  motion 
with  radiation  pressure  acceleration  in  Earth-centered  inertial 
coordinates  (ECI)  are  given  by 


[ax] 


0  -a3  a2 

a3  0  -a  1 
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is  the  skew-symmetric  matrix  representation  of  the  cross 
product  for  any  general  3x1  vector  a.  The  quaternion 
kinematics  equation  is  given  by 


q  =  t-H(q)ro  (17) 

where  co  is  the  component’s  angular  velocity.  The  angular 
velocity  dynamic  equation  can  be  written  as 


(b  =  J_1  (m  -  [co  xJ/(o)  (18) 

where  J  the  inertia  tensor  for  the  SO  and  M  are  any  external 
applied  torques.  The  radiation  pressure  moments  can  be 
calculated  by  considering  that  the  forces  act  through  the  center 
of  each  facet 


-r7+a7 


(12) 


where  the  terms  fi  represents  the  gravitational  parameter  of  the 
Earth,  a^2  is  the  gravitational  perturbation  due  to  non- 
symmetric  distribution  of  mass  along  the  lines  of  latitude,  and 

^ facets 

2l ^1  -  ^a7  represents  the  acceleration  perturbation  due  the 

i= 1 

various  radiation  pressures  and  summed  over  all  the  surfaces. 
Details  regarding  the  calculation  of  this  last  term  can  be  found 
in  [3]. 

The  attitude  matrix  for  each  component  of  the  SO  can  be 
written  as  a  function  of  the  component’s  quaternion  by 


A  =  S{qJ^{q)  (13) 


where 


H(q)  = 


C?473x3  +  P5  xt 
~\T 


(14) 


^(q)= 


-%T 


(15) 


M  =  m  Z[r'S  xf(q)a^  (19) 

M 

where  is  the  location  of  the  geometric  center  of  each  facet 

with  respect  to  the  center  of  mass  of  the  SO  in  body 
coordinates  and  A( q)  is  the  attitude  matrix  calculated  by  the 
quaternion  q. 

The  surface  properties  are  assumed  to  be  constants,  and  so 
their  dynamics  equation  can  be  written  as 


p  =  0  (20) 

IV.  Unscented  Kalman  Filter 

In  this  paper,  the  quaternion-based  unscented  Kalman  filter 
(UKF)  of  Crassidis  and  Markley  [6,7]  is  employed  with  the 
attitude  state  errors  represented  as  error  Generalized  Rodrigues 
Parameters  (GRPs)  [8].  The  augmented  state  is  the  SO’s 
orientation  and  rotation  rate,  and  the  SO’s  position,  velocity, 
and  surface  parameters  given  by  (1 1). 

The  parameter  propagation  equations  in  (12),  (17),  (18), 
and  (20)  can  be  written  in  the  general  state  function  which 
gives  the  deterministic  part  of  the  stochastic  model 

(2i) 


and 


where  wk  is  the  process  noise  vector 

In  this  approach,  the  observation  vector  includes  the 
apparent  magnitude  as  computed  with  (10),  as  well  as 


astrometric  measurements  of  the  right  ascension  and 
declination  of  the  SO: 

y  =  \m  object  RA  Dec]  (22> 

where  y  is  a  1^3  vector  of  real  numbers.  The  general 
measurement  function  used  in  the  estimation  filter  is: 

yk=h{x~k,vk)  (23) 

where  vk  is  the  measurement  noise  vector. 

V.  Simulations 

The  particular  example  will  be  that  of  simultaneously 
estimating  the  position,  velocity,  attitude,  angular  rates,  and 
surface  parameters  of  a  HAMR  object  in  geosynchronous  Earth 
orbit  (GEO).  Table  I  lists  the  initial  truth  state,  the  initial 
estimated  state,  and  the  initial  uncertainty. 

The  shape  is  defined  as  a  cube  with  1  -m  sides  and  a  mass  of 
2  kg.  Each  surface  of  the  cube  is  coated  with  the  same  BRDF 
surface  parameters.  The  orbit  was  set  to  geosynchronous  (< a  = 
42364.16932  km,  e  «  0,  /  *  30°,  M0  =  91°,  on  =  0,  Q  =  0). 
Observations  were  simulated  starting  at  2010  Mar  15  at 
4:00:00  UT,  1800  observations  every  2  s  (for  a  total  of  1  hour), 
with  an  observation  site  corresponding  to  the  top  of  Haleakala 
on  Maui  (latitude  =  20.71  deg,  longitude  =  -156.26  deg, 
altitude  =  3.0586  km).  The  Thermal  Radiation  Pressure  (TRP) 
parameters  of  each  surface  were  C  =  9000  J/K,  K  =  25.5  W/K, 
and  Tbody  =  243.5  K.  The  measurement  noise  is  0.1  mag  for  the 
brightness  observation  and  10  arc-sec  in  the  right  ascension 
and  declination  observations. 


TABLE  I.  Tests  #l/#2  Augmented  State  Setup 


Value  in 
State 

Initial  Truth 

Initial  Estimate 

Uncertainty 

0.754 

0.695 

q 

0.133 

0.000 

0.134 

0.010 

3.33  deg 

0.643 

0.706 

0.00200 

0.00212 

w  (rad/s) 

-0.00100 

-0.00106 

1.16  x  10*4 

0.00500 

0.00506 

-739.4 

-789.4 

r1  (km) 

36682.9 

36732.9 

100 

21178.9 

21278.9 

-3.0669 

-3.0169 

v7  (km/s) 

-0.0464 

0.0536 

0.10 

-0.0268 

-0.0768 

n 

150 

140 

10 

P 

0.40 

0.30 

0.10 

d 

0.70 

0.80 

0.10 

In  the  first  two  tests,  the  physically-consistent  BRDF/Solar 
Radiation  Pressure  (SRP)/TRP  model  was  used  to  generate  the 
truth.  In  the  first  test,  this  same  model  was  used  in  the  state 
function  of  the  UKF  estimation  while  in  the  second  test,  a 
simplified  SRP  model  (where  the  BRDF  has  a  Lambertian 
diffuse  component  and  a  mirror-like  specular  component)  and 
no  TRP  model  was  used  in  the  state  function.  Figs.  2-4  plot  the 
measurements  over  the  observation  period. 

Figs.  5-7  plot  selected  components  of  the  estimated  state 
and  covariance  as  a  function  of  time  for  Test  #1  while  Figs.  8- 
10  plot  the  same  for  Test  #2. 


Fig.  2.  Brightness  in  magnitudes  as  function  of  time 


Fig.  3.  Right  Ascension  as  function  of  time 


Fig.  4.  Declination  as  function  of  time 


10 


10 


Fig.  5.  Attitude  difference  from  truth  and  3 -a  error  bounds  for  Test  #1 


Fig.  8.  Attitude  difference  from  truth  and  3-a  error  bounds  for  Test  #2 


Fig.  6.  Angular  rate  difference  from  truth  and  3-o  error  bounds  for  Test  #1 


Fig.  9.  Angular  rate  difference  from  truth  and  3-o  error  bounds  for  Test  #2 


Time  (min) 


i - 

1 - 

1 - 

- / 

r”  '  . . 

i _ 

i _ 

i _ 

i _ 

0  10  20  30  40  50  60 

Time  (min) 


Fig.  7.  Surface  parameter  difference  from  truth  and  3-a  error  bounds  Fig.  10.  Surface  parameter  difference  from  truth  and  3-a  error  bounds 

for  Test  #  1  for  T est  #2 


TABLE  II. 


Tests  #l/#2  Results 


Value  in 

Test  #1 

Test  #2 

State 

Difference 
from  Truth 

Uncertainty 

Difference 
from  Truth 

Uncertainty 

0.152 

0.131 

1.628 

0.126 

angle 

-0.267 

0.106 

-1.365 

0.110 

(deg) 

0.828 

0.391 

4.898 

0.376 

CO 

-0.273 

1.467 

-3.288 

1.762 

(deg/hr) 

3.253 

3.131 

27.184 

2.972 

-0.158 

0.509 

-6.898 

0.547 

29.4 

22.4 

28.1 

22.4 

r'(km) 

-75.5 

57.4 

-12  A 

57.4 

-46.6 

35.5 

-44.5 

35.5 

0.00694 

0.00524 

0.00663 

0.00524 

v7  (km/s) 

-0.00046 

0.00043 

-0.00045 

0.00043 

-0.00033 

0.00026 

-0.00032 

0.00027 

n 

0.775 

1.493 

2.796 

1.462 

P 

-0.0053 

0.0032 

-0.0164 

0.0035 

d 

-0.0015 

0.0031 

0.0076 

0.0035 

Table  II  lists  the  difference  between  the  final  estimated  state 
and  truth,  and  the  1-g  uncertainty  as  calculated  by  the  final 
covariance  matrix  for  both  tests.  The  attitude  quaternion 
difference  has  been  converted  to  the  equivalent  roll/pitch/yaw 
Euler  angle  differences. 

Of  particular  note  in  Test  #1  is  that  all  parameters  in  the 
state  are  observable  by  the  filter  and  converge  to  the  correct 
values  (i.e.  all  parameters  to  within  3-g  of  the  quoted 
uncertainty).  In  Test  #2,  however,  due  to  the  dynamics  model 
mismatch,  some  parameters  in  the  state  are  not  converging,  and 
the  differences  from  truth  of  the  final  state  for  these  parameters 
are  well  outside  3-g  of  the  quoted  uncertainty  (e.g.  attitude 
Euler  angle  components  >  10-g  from  truth).  Some  of  the 
parameters  are  unaffected  (e.g.  position  and  velocity)  by  the 
model  mismatch  over  the  one  hour  time  scale  sampled,  but 
would  undoubtedly  be  affected  over  a  longer  time  span. 

In  the  next  two  tests,  the  same  initial  conditions  with  the 
physically-consistent  BRDF/SRP/TRP  used  to  both  generate 
the  truth  and  in  the  estimation  filter  as  in  Test  #1  are  duplicated 
except  that  the  uncertainties  on  the  surface  parameters  are 
increased  by  a  factor  of  3  from  those  shown  in  Table  I. 
Additionally,  in  Test  #4,  the  difference  from  the  truth  and  the 
uncertainties  of  the  orbital  and  attitudinal  state  parameters  are 
decreased  by  a  factor  of  10. 

Figs.  11-15  display  how  the  increased  initial  covariance  for 
the  surface  parameters  results  in  filter  divergence  by  plotting 
all  the  state  parameter  differences  from  truth  and  the  3-g 
uncertainty  as  a  function  of  time  for  Test  #3.  In  this  test,  even 
the  position  and  velocity  begin  to  diverge.  In  contrast,  when 
the  orbital  and  attitudinal  state  parameter  uncertainties  are 
decreased  and  the  estimation  rerun,  as  in  Test  #4,  the  filter  is 
able  to  converge  (although  outside  the  3-g  bound),  despite  the 
increased  surface  parameter  uncertainty  as  shown  in  Fig.  16. 


Fig.  1 1 .  Attitude  difference  from  truth  and  3 -a  error  bounds  for  Test  #3 


Fig.  12.  Angular  rates  difference  from  truth  and  3 -a  error  bounds  for  Test  #3 


Time  (min) 

Fig.  13.  Position  difference  from  truth  and  3 -a  error  bounds  for  Test  #3 


TABLE  III.  Tests  #3/#4  Results 


Time  (min) 


Fig.  14.  Velocity  difference  from  truth  and  3-o  error  bounds  for  Test  #3 


Time  (min) 


Fig.  15.  Surface  parameter  difference  from  truth  and  3-o  error  bounds 
for  Test  #3 


Fig.  16.  Surface  parameter  difference  from  truth  and  3 -a  error  bounds 
for  Test  #4 


Value  in 

Test  #3 

Test  #4 

State 

Difference 
from  Truth 

Uncertainty 

Difference 
from  Truth 

Uncertainty 

Euler 

40.730 

0.119 

-0.457 

0.119 

angle 

-10.058 

0.130 

-0.357 

0.100 

(deg) 

-14.387 

0.194 

0.728 

0.248 

36.226 

0.592 

-1.506 

0.918 

CO 

-112.338 

1.087 

3.807 

1.599 

(deg/hr) 

434.556 

0.545 

-1.151 

0.413 

-279.5 

21.5 

-3.5 

9.4 

r7  (km) 

715.9 

54.9 

8.9 

24.0 

442.6 

34.0 

5.5 

14.8 

-0.06546 

0.00503 

-0.00078 

0.00220 

v7  (km/s) 

-0.00033 

0.00043 

-0.00035 

0.00043 

-0.00040 

0.00027 

-0.00023 

0.00027 

n 

-208.2 

4.738 

-2.666 

1.630 

P 

-0.4960 

0.0027 

0.0248 

0.0017 

d 

0.2388 

0.0010 

-0.0243 

0.0023 

Table  III  lists  the  difference  between  the  final  estimated 
state  and  truth,  and  the  1-g  uncertainty  as  calculated  by  the 
final  covariance  matrix  for  Test  #3  and  Test  #4.  The  attitude 
quaternion  difference  has  again  been  converted  to  the 
equivalent  roll/pitch/yaw  Euler  angle  differences.  In  Test  #3, 
the  filter  is  suffering  from  information  dilution.  There  are 
simply  too  many  parameters  in  the  state  whose  uncertainties 
are  large.  When,  however,  as  in  Test  #4  the  uncertainties  for 
only  a  limited  number  of  state  parameters  are  large,  the  filter  is 
able  to  converge  to  close  to  the  truth. 

VI.  Conclusions 

Physically  consistent  BRDF,  SRP,  and  TRP  modeling 
enables  astrometric  and  photometric  data  fusion  for  the 
purposes  of  simultaneously  estimating  orbital,  attitudinal,  and 
surface  parameters  of  a  space  object.  Other  parameters 
associated  with  the  space  object,  such  as  the  mass  and  shape, 
could  also  be  added  to  the  augmented  state  vector  and 
estimated  provided  that  their  effect  on  the  system  is  observable 
with  an  appropriate  measurement  model.  Caution  is  warranted, 
however,  concerning  model  mismatches  (both  in  the  dynamics 
and  the  BRDF  model)  and  the  possibility  of  information 
dilution  when  too  many  quantities  with  large  uncertainties  are 
estimated  at  once. 
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